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Abstract 

Within the framework of density functional theory (DFT), the total energy of crys¬ 
tal structures is calculated at zero temperature. Herein, we briefly discuss the DFT-based 
lattice-dynamics approach for computing crystal free energy, the quantity needed in various 
non-zero-temperature contexts. We illustrate this well-established approach by examining 
the temperature-dependent thermodynamic stability of several crystalline materials, in¬ 
cluding ZrOa, HfOa, KBH 4 , and Zn(BH 4 ) 2 . 


1 Introduction 

Density functional theory (DFT) QH is a powerful numerical method in various disciplinaries, 
e.g., computational chemistry and computational materials science. In a typical calculation 
with DFT, the total energy of a static atomic (or more precisely, ionic) structure {r*} — here 
Vi is the coordinate of the ion — of a solid structure of N atoms, is evaluated. The forces 
exerting on these ions can also be computed with the Hellmann-Feynman theorem and then 
be used to relax the structure, i.e., to locate a local minimum of the potential energy surface 
in the conhgurational space. Some materials properties, e.g., band gap, can then be calculated 
on top of the relaxed structures obtained. Because temperature is not included in the DFT 
formalism, all of these properties, when calculated with DFT, are valid at zero temperature 
(T = 0). Of particular interest are, however, certain materials properties at some sort of hnite 
temperatures, i.e., T > 0. Technically, it is possible to evaluate some properties at non-zero 
temperatures, of course, with certain computational overhead. 

This contribution describes a method, which is based on lattice dynamics [^, for estimating 
the free energy of crystals at non-zero temperatures. Presumably, this is one of the most 
desirable hnite-temperature quantities as it determines pretty much the essential characters 
of a crystal [5||^. To be more precise, the free energy difference between systems and the 
gradients of the free energy with respect to thermodynamic variables are of practical interest. 
For instance, among various possible crystal structures, that with lowest-possible free energy 
is considered to be thermodynamically stable. For a chemical reaction, the free energy change 
is an useful parameter, indicating whether the reaction is spontaneous or not, and if not, how 
much energy is needed to make it occur. 

The Gibbs free energy of an equilibrium (relaxed) structure {r?} at constant number of 
particles is dehned as G{P,V,T) = U — TS -\- PV where U is the internal energy, S is the 
entropy, P is the pressure, and V is the volume (per unit cell, for examples) of the crystals. 
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Strictly, the ions are not frozen but perform certain vibrations in the neighborhood of {r^}. For 
this reason, the internal energy may be expressed asU = 17o + hvib where Uq is the static internal 
energy at {r?} and ?7vib is the vibrational internal energy [^. Technically, Uq is the total DFT 
energy of the relaxed structure {r^} while PV can also readily be computed. The treatment 
for the remaining terms, i.e., the vibrational free energy dehned as F’vib(T') = U^ib — TS, which 
is considered in details herein, is based on the calculations of the vibrational frequency spectra 


showing some advantages and shortcomings which will also be discussed in this contribution. 


of the examined solids. This approach is well-established and has widely been used [8-16 


2 Vibrational free energy 

To focus on FVib(r), we put the PV term of G aside, and consider the Helmholtz free energy 
F{T) = U — TS = Uo + Fvib(T). The formal dehnition of F{T) is 

F{T) = -kBT In Z (1) 


where k-Q is the Boltzmann’s constant, and Z is the canonical partition function dehned by 


Z= djrjexp 


ksT 


( 2 ) 


Here, the integral is taken over all the possible atomic conhguration {rj} of the system, and 
H = T + ^ where T is the kinetic energy of the ionic lattice and $ = <h({rj}) is the potential 
of the ion structure. In principles, the very-high-dimensional integral appearing in Eq. Q can 
not be performed directly. It may, however, be sampled from a thermodynamic ensemble by 
a number of techniques, many of them are based on Monte Carlo like the Bennett acceptance 
ratio method [l7| . For a better summary of such methods, readers are referred to Ref. [l^ . 

The lattice-dynamics approach for F^i\,{T) is started from the harmonic approximation of 
the potential energy <I>({rj}), i.e., 
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$({ra) = $o({r°}) + E 




2 i^,fs dri,adr,^p 




(3) 


In this expansion, Ui^a is the displacement of ion i (j) along the a {(5) direction. It 

should be noted here that $0 = the internal energy of the frozen ionic structure {r?} and 
t^vib = T + ^ (3), the displacements are assumed to be small, 

and hence unharmonic terms like those arisen from substructure rorations and the translations 
(these motions may be possible in some highly dynamical solids like complex borohydrides 
|19f|21] ) are not included. In the quantum description of the lattice vibrations, a set of normal 
coordinates {qu} and pi^ = —ifid/dqy {fi is the Planck constant) are introduced so that 


H = Uq + T + 
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where {cOi,} are the vibrational frequencies. The eigenenergies of H are then 


3N 


= Uq + ^ [uu + - \ huj, 


u=l 


(5) 
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in which n = {ni, n 2 ,nsAr}. By inserting en into Q, one gets Z = exp(—C/o/^BT")^vib with 
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Eq. Q implies that Fvib(r) is determined by Fvib(T) = —/sBTlnZvib, or 
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( 6 ) 


(7) 


At the limit of T —)• 0, -Fvib(7~') approaches the zero-point energy, which is given by 


3N 


Ezp = 




U=1 


( 8 ) 


In case the spectrnm of uj is continnons, Evib(iE) is written in terms of the normalized density 
of phonon states g{u}) as 


POO 

Evih{T) = SNk^T / dujg{u})ln 

Jo 


2 sinh 
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(9) 


Thns, given that the density of phonon states g{uj) is determined by whatever method, the 
vibrational free energy F^i\^(T) can readily be calcnlated. 


3 Calculations of g{uj) with DFT 


The density of phonon states g{oj) can be calcnlated nnmerically by some methods. Here, we 
examine the prescriptions made available within the DFT-based framework, leaving aside other 
methods, e.g., that starts from the velocity antocorrelation fnnction which can be accessible 
with molecnlar dynamics. 

Calcnlations of g{uj) with DFT involve the determination of phonon dispersion relation 
— here q is the phonon wave vector — by diagonalizing the dynamical matrix, dehned 


as 


7^0/3 (D;q) 
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( 10 ) 


where m* and rrij are the mass of ions i and j, respectively. The snm in (10) is taken over the 
index I of the translational images of the crystal cell dne to the periodicity. With the dynamical 
matrix D introdnced, the eqnation of motion of the vibrations is written as 


'^Dafj{ij;q)e0{j,ciu) = mja;^(q)ea(i, qi^) (11) 

where ea{i, q^^) are the coefficients of the transformation from Ui^a to j^. Thns, the problem 
of calcnlating the phonon dispersion ujy{(i) of the normal phonon mode is rednced to the 
eigenvalne problem of Zla^(ij;q), which can be solved nnmerically. 

In practice, DQ,^(ij;q) and Wi/(q) are compnted by either the frozen-phonon method (also 
referred to as the direct approach) |22] or the linear-response method |23f|^ . The former 
relies on evalnating the interatomic force constants of a large snpercell in which the ions are 
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perturbed from their equilibrium positions. The size of the supercell may be large, depending 
on the commensurability of the perturbations with the equilibrium cell. Then, the obtained 
matrix of force constants is reduced to the dynamical matrix at certain values of q which is 
then be diagonalized to obtain 0 Ju{q). In this method, a DFT package like VASP 26-29 or 
SIESTA is needed as a force calculator, while an additional software, e.g., PHONON 22 , 
PHONOPY |31|, and PHON |32| , is used for the pre- and post-processing calculations. 

The linear-response method for calculating oj^{c^ and g{uj) is based on the perturbation 
theory within DFT |23] . In this method, the dynamical matrix is evaluated from the knowledge 
of the hrst-order response of the wave function with respect to the lattice potential perturbations 
|24[|25| . Different from the frozen-phonon method which needs supercell to be constructed, the 
dynamical matrix is separately evaluated by linear-response calculations on the unit cell at 
various q points. Available packages that support the linear-response method for calculations 
of a;jy(q) are, but not limited to, ABINIT and QUANTUM ESPRESSO p4] . 

Given that the phonon dispersion is determined, the density of phonon states g{oj) 

and the vibrational free energy can be straightforwardly calculated, as supported by 

all of the codes/tools mentioned above. In addition, calculated Wi./(q) is also very useful to 
examine the dynamical stability of crystal structures, especially those predicted theoretically 
T 0 H 13 15,35 . In particular, if a structure corresponds to a saddle point of the potential energy 


surface, there are at least one phonon modes with imaginary frequencies, and one can follow 
these modes to end up at one or more lower-energy, dynamically stable structures 
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4 Free energy and thermodynamic stability 

In this Section, we examine two examples of using the calculated free energy for some discussions 
of the thermodynamic stability of crystals. 


4.1 Temperature-driven structural phase transition 


At a given temperature/pressure condition, the lowest Gibbs free energy structure of a certain 
crystal is considered to be thermodynamically stable. Gonsequently, the structural phase tran¬ 
sitions of a crystal may be traced by examining the Gibbs free energy calculated for its possible 
13-^. Here, we limit at P = 0, thus the Helmholtz free energy F{T) will be used 


structures 

to estimate the critical temperature Tq of the phase transitions of some representative crystals. 

Three crystalline solids which are used for this demonstration, are summarized in Table 
Both zirconia Zr02 and hafnia Hf02 are known |36p8 -40] to transform from a monoclinic P2i/c 
structure at low temperatures to a tetragonal PA 2 /nmc structure at higher temperatures. These 
phase transitions occur at Tq ~ 1500K and Tq ~ 2200K for zirconia and hafnia, respectively. 


Table 1: Low- and high-temperature phases of some crystalline materials examined in this 
work. _ 


Materials 

Low-T structure 

High-T structure 

Tc Ref. 

Zr 02 

P 2 i/c 

P4:2/nmc 

1500K 

36 


Hf 02 

P2ilc 

PA2/nmc 

2200K 

36 


KBH 4 

PA2/nmc 

Fm3m 

197K 

37 
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Figure 1: Helmholtz free energies calculated for the high-T phases of Z 1 O 2 , Hf 02 and KBH 4 
with respect to their low-T phases. Details on these phases can be found in Table 


Likewise, potassium borohydride KBH 4 crystallizes in a tetragonal P 42 /nmc structure below 
Tc — 197 K while adopting a cubic Fm3m structure above this point [^. The Helmholtz free 
energy F{T) has been calculated for all of these structural phases, using VASP and PHONOPY. 
In Fig. [^we show the free energy of the high-temperature (high-T) structure of each crystal 
with respect to that of the low-temperature (low-T) structure. Clearly, the low-T and high- 
T structures of zirconia, hafnia, and potassium borohydride are correctly described by the 
calculations. The phase transitions of zirconia and hafnia are predicted to be 1600K and 2100 
K (T^ , which are in relatively good agreements with the known facts. On the other hand, the 
critical temperature predicted for potassium borohydride is about 550 K, which is considerably 
higher than the correct value (a similar observation was also reported in Ref. |14]). 

The discrepancy observed for Tc of potassium borohydride may be understandable, given 
that the lattice-dynamics approach for calculating Tvib(T) assumes a series of approximations 
and simplifications. Moreover, as mentioned above, complex borohydrides like KBH 4 are highly 
dynamical, i.e., the BH 4 groups may jump and/or rotate 19 -^. Such the motions can not be 


captured by the harmonic approximation, which assumes small displacements. Next, the lattice- 
dynamics approach, as described in this contribution, is basically an extrapolation procedure, 
i.e., Uq and g{u}) of the high-T cubic Fm3m phase, which exists at T ~ 197K, are calculated 
at OK and then, the relevant energetic quantity, i.e., T(T), is extrapolated back to non-zero 
temperatures using Eq. (§. Finally, the energy difference between the high-T and low-T phases 
of potassium borohydride is very small (~ 3 meV/atom, more than one order smaller than those 
of zirconia and hafnia), suggesting that this numerical scheme may still not be adequate for a 
marginal case like KBH 4 . 
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Figure 2: Free energy change AF[T) of reaction (12). 


4.2 Thermodynamic stability w.r.t. formation/decomposition 


For another demonstration, we consider a chemical reaction of which the reactants and products 
are all solids |41[|42| 

2 NaBH 4 + ZnCh Zn(BH 4)2 + 2NaCl. (12) 

This is the synthesis route of Zn(BH 4 ) 2 , as reported in Refs. |41}|42] . However, the stability of 
Zn(BH 4)2 is still in a debate, and a the presence of pure Zn(BH 4)2 is not conclusive. While the 
structure of the synthesized Zn(BH 4)2 samples has yet been resolved, a number of low-energy 
structures have been proposed via DFT calculations flO 43, 4^ . As one of the simple tests 


that may be performed on these hypothetical structures, the free energy change AF{T) of the 


reaction (12) can be calculated. A positive value of AF{T) suggests that the examined structure 
of Zn(BH 4)2 is unstable, while a negative value of AF(T) may be a supporting evidence for 
Zn(BH4)2. 

We show in Fig. (|^ the free energy change AF (T) calculated for the reaction (12), using the 
currently known most-stable structures of the crystals involving. They are the /4i22 structure 
of Zn(BH 4)2 10 , the P42/nmc structure of NaBH 4 37 , the Fm3m structure of NaCl, and the 


, the P4:2/nmc structure of NaBH 4 

Pna2i structure of <5—ZnCb |^. A conclusion which may be drawn from Fig. ([^ is that up 
to 500 K, the formation of Zn(BH 4)2 according to the reaction (12) is energetically permissible. 
It is however worth noting that the thermodynamic stability of Zn(BH 4)2 should be examined 
with respect to all the possible formation/decomposition chemical reactions, the practice which 
is well beyond the scope of this contribution. 


5 Conclusions 

In summary, Gibbs and Helmholtz free energies of crystals may be estimated by lattice dy¬ 
namics within the harmonic approximation at the level of density functional theory. This 
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computational scheme produces reasonable results for many crystals but there are also other 
solids with which this approach should be used with great cautions. Although some treatments 
for the anharmonic terms of the free energies are also currently available, they should also be 
material-dependent in the similar fashion with this scheme. 
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